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1  Introduction 


Background 

Shocks  in  fluids  result  from  fluid  flow  that  is  more  rapid  than  the  ^ecd  of 
a  compression  wave.  TTius  there  is  no  means  for  the  flow  to  adjust  gradually. 
Pressure,  velocity,  and  temperatures  change  abruptly,  causing  severe  fatigue 
and  component  destruction  in  military  aircraft  and  engine  turbines.  This 
problem  is  not  limited  to  supersonic  aircraft;  many  parts  of  subsonic  craft  are 
supersonic.  For  example,  the  rotors  of  helicopters  have  supersonic  regions  as 
do  the  blades  of  the  turbine  engines  used  on  many  crafts.  The  shock  is  formed 
as  the  flow  passes  from  supersonic  to  subsonic  or,  in  the  case  of  an  oblique 
shock,  as  the  result  of  a  geometric  transition  in  supersonic  flow.  Wind  tunnels 
are  limited  in  the  Mach  numbers  they  can  achieve  and  testing  is  expensive; 
thus  design  relies  upon  numerical  modeling.  In  hydraulics  the  equivalent 
shocks  are  referred  to  as  hydraulic  jumps,  surges,  and  bores.  Here,  for 
example,  it  is  important  to  determine  the  ultimate  height  of  water  resulting 
from  a  dam  break  or  the  insertion  of  a  bridge  in  a  fast- flowing  river. 

The  compressible  Euler  equations  describe  these  flow  fields  and  are  solved 
numerically  in  discrete  models.  These  partial  differential  equations  implicitly 
assume  a  certain  degree  of  smoothness  in  the  solution.  Models,  therefore, 
have  great  difficulty  handling  shocks.  One  method  to  avoid  solving 
numerically  across  the  shock  is  to  track  the  shock  and  impose  an  internal 
boundary  there.  This  method  is  called  "shock-tracking."  On  the  other  hand 
the  sharp  resolution  of  the  shock  can  be  forfeited  and  allow  for  0(1)  error  at 
the  transition.  This  is  referred  to  as  "shock-capturing,"  as  originated  by 
von  Neumann  and  Richtmyer  (1950),  and  is  now  the  most  common  technique 
used  in  engineering  practice. 

Great  care  must  be  undertaken  to  make  sure  the  errors  are  local  to  the 
shock.  Otherwise  the  shock  location  and  speed  will  be  incorrect.  It  is 
important  that  the  discrete  numerical  operations  preserve  the  Rankinc-Hugoniot 
condition  (Anderson,  Tannehill,  and  Fletcher  1984)  resulting  from  integration 
by  parts.  While  this  should  result  in  reasonable  shock  speed  and  location, 
discrete  models  commonly  suffer  from  numerical  oscillations  near  the  shock. 

There  are  many  proposed  methods  used  to  reduce  these  oscillations.  The 
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basic  theme  is  to  cleverly  apply  artificial  diffusion  as  a  result  of  flow 
parameters.  Many  of  these  methods  do  not  preserve  the  original  equations 
within  the  shock  due  to  this  added  diffusion.  Hughes  and  Brooks  (1982)  have 
approached  this  problem  within  the  finite  elements  method  by  the  development 
of  a  single  test  function  that  reflects  the  speed  of  fluid  transport  (the  SUPG 
scheme,  Streamline  Upwind  Petrov-Galerkin)  to  be  applied  to  the  entire  partial 
differential  equation  set.  In  this  manner  the  model  is  consistent  even  at 
discrete  scales.  Its  application,  thus  far,  has  been  only  to  the  very  simple  case 
of  Burgers*  equations. 

In  this  report  a  two-dimensional  (2-D)  finite  element  model  for  the  shallow- 
water  equations  is  produced  using  an  extension  of  the  SUPG  concept,  but  rely¬ 
ing  upon  the  characteristics  of  the  advection  matrix  (transport  as  well  as  the 
free-surface  wave  speeds)  to  develop  the  test  function  to  be  applied  to  the 
coupled  set  of  equations.  The  shallow-water  equations  are  a  direct  analogy  to 
the  Euler  equations  with  the  depth  substituted  for  density  and  dropping  the 
energy  equation.  This  equation  set  is  ideal  for  testing  numerical  schemes  for 
the  Euler  equations.  The  model  developed  can  reproduce  supercritical  and 
subcritical  flow  and  is  shown  to  reproduce  very  difficult  conditions  of 
supercritical  channel  transitions  and  preserve  the  Rankine-Hugoniot  conditions. 

For  simple  geometries,  analytic  and  flume  results  are  compared  with 
approaches  for  shock-capturing  and  shock  detection.  A  trigger  mechanism  that 
turns  on  the  capture  schemes  in  the  vicinity  of  shocks  and  the  characteristic 
upstream  weighted  test  function  are  tested. 

Tue  results  of  this  research  are  an  algorithm  and  program  to  represent 
hydraulic  jumps  and  oblique  jumps  in  2-D  for  shallow  flow.  The  code, 
HIVEL2D,  is  a  general-purpose  tool  that  is  applicable  to  many  problems  faced 
in  high-velocity  hydraulic  channels,  notably,  in  the  calculation  of  the  ultimate 
water  surface  height  around  bridges,  channel  bends,  and  confluences  subjected 
to  supercritical  flow  or  due  to  surges  caused  by  sudden  releases  or  dam  failure. 
The  algorithm  itself  is  applicable  outside  the  field  of  hydraulics  as  well  to 
complex  aerodynamic  flow  fields  containing  shocks. 


Basic  Equations 

The  basic  equations  that  are  addressed  are  the  2-D  shallow-water  equations 
given  as: 


dt  dx  dy 


+  //  =  0 


(1) 
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where 


where 

t  =  time 

x,y  =  Horizontal  Cartesian  coordinites 
h  =  depth 

p  =  x-discharge  per  unit  width,  uh 
q  =  y-discharge  per  unit  width,  v/i 
g  =  acceleration  due  to  gravity 
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(da  dv\ 

ai2  =  02i  =  pv'-+-J 

-  3v 
=  2pv  ^ 

p  =  fluid  density 

V  =  kinematic  viscosity  (turbulent  and  molecular) 
w,v  =  velocity  in  x,  y  directions 
dz 

z  =  bed  elevation 
n  =  Manning’s  coefficient 
c  =  1.0  metric,  1.486  non-SI 
dz 

^2“3y 

These  equations  neglect  the  Coriolis  effect  and  assume  the  pressure  distribution 
is  hydrostatic,  and  the  bed  slope  is  assumed  to  be  geometrically  mild  though  it 
may  be  hydraulically  steep.  These  equations  apply  throughout  the  domain  in 
which  the  solution  is  sufficiently  smooth.  Now  consider  the  case  for  which 
the  solution  is  not  smooth. 


Shock  equations 

In  this  section  we  first  derive  the  jump  conditions  in  one  dimension  (l-D) 
with  no  dissipative  terms,  i.e.,  friction  or  viscosity.  We  show  that  as  a  result 
of  the  discontinuity  of  the  jump,  the  shallow-water  equations  should  conserve 
mass  and  momentum  balance  but  will  lose  energy.  Furthermore,  if  there  is  no 
discontinuity,  energy  too  will  be  conserved.  Later  the  jump  relations  are 
extended  to  2-D. 

This  derivation  relies  upon  the  work  of  Stoker  (1957)  and  Keulegan  (1950) 
following  a  fluid  element  through  a  moving  jump.  Figure  1  defines  these 
features. 

If  our  coordinate  system  is  chosen  to  move  with  the  jump  at  speed  V^,  then 
we  may  use  the  following  term  definitions. 
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Figure  1 .  Definition  of  terms  for  a  moving  hydraulic  jump 


Uq  -  velocity  at  section  0 
iiQ  =  depth  at  section  0 
=  velocity  at  section  1 
hj  -  depth  at  section  1 

Vq  =  Uq'  velocity  at  section  0  relative  to  the  jump 

Vy  =  relative  velocity  at  section  1 

Mass.  Now  following  an  infinitesimal  fluid  clement  in  our  moving  coordi¬ 
nate  system  we  know  that  mass  is  conserved  so  we  have, 

(2) 

dt 


where  p,  the  fluid  density,  is  a  constant  here. 
Across  the  element  we  have 
P(V'i/«i  -  Vq/'o)  =  0 
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or 


=  ViQ  =  q 


(4) 


where  q  is  the  relative  discharge.  Equation  3  may  be  written  in  a  fixed  grid  as 
V^lh^  =  lUh]  (5) 

where  the  symbol  [  J  implies  the  jump  in  the  quantities  across  the  discon- 
tinuuy,  e.g.,  Vti]"  hj  -  Hq. 

Momentum.  In  the  same  manner  we  show  that  momentum  is  conserved 
by: 

5^  =  p(Po  -  Pi)  (6) 

at 

P  -  (VQhQ)UQ]  =  p  (Pq  -  ^i)  ('7) 

where 

D  1/2 

D  1  /.2 


^  (t/i  -  I/o)  "  Po  -  e,  (8) 

v/hich  in  a  fixed  coordinate  system  would  be: 

[Uh]  =  [f/2  /,  +  PJ  (9) 

Energy.  Now  consider  the  case  of  mechanical  energy  E  as  it  passes 
through  the  discontinuity 

■1  ''I'M''?  -  1  Vot'i?  * 

"  ^0^0  *  ^1^1  *  ^o^oj  (’.0) 

If  we  multiply  (8)  by  add  this  to  (10)  using  the  relationship  for  q  we  have: 
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j  qiV^  -V^-2  -  P^^Vq) 


=  P 


~  qiVl  -  Vq)  (Vi  -  Vq)  -  2  (PiVj  -  PqVq) 


(U) 


Now  substituting  Equation  8  yields 


dE 

~di 


P 


\  (Pq  -  ^l)  (Vj 


-Vq)*2  (PjV, 


-  Po'^o) 


(12) 


or,  finally 


dE  _  pqg  (^'o  -  ^'l) 
dt  4  /jq/i  j 


(13) 


so  that  the  shallow-water  equations  lose  energy  at  the  jump,  and  it  is  propor¬ 
tional  to  the  depth  differential  cubed.  If  the  depth  is  continuous,  no  energy 
will  be  lost. 

While  mathematically  an  energy  gain,  dEldl  >  0,  through  the  jump  is  a 
possibility,  physically  it  is  not.  If  we  restrict  ourselves  to  energy  losses 
through  the  jump,  there  are  two  possible  cases. 

a.  Case  1:  Vq,V^<0  implies 

hQ>  hj.  The  jump  progresses 
downstream  through  our  fluid 
clement  (Figure  2). 

b.  Case  2:  Vq,  Vj  >  0  implies 
Hq  <  h j.  The  fluid  passes 
downstream  through  the  jump 
(Figure  3). 

Here  we  have  arbitrarily  chosen 
the  flow  to  be  from  left  to  right,  if 
we  had  chosen  the  opposite  direction 
one  would  simply  have  horizontal 
mirror  images  of  Figures  2  and  3.  A  fluid  particle  that  is  about  to  be  swept 
into  or  caught  by  the  jump  is  considered  in  "front"  of  the  jump.  A  fluid 
element  that  has  passed  through  the  jump  is  now  "behind"  it.  Therefore,  we 
may  conclude  that  the  water  level  is  lower  in  front  of  the  jump  than  it  is 
behind  the  jump. 

In  order  to  calculate  the  wave  speed,  it  is  convenient  to  chexise  Uj  =  0. 


Figure  2.  Example  of  Case  1;  jump 

passes  downstream  through 
the  fluid  element 
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— > 

Front 

Back 

Figure  3.  Example  of  Case  2;  the  fluid 
passes  downstream  through 
the  jump 


This  is  completely  arbitrary,  and  this 
form  also  produces  an  easy  test  case 
that  will  be  eventually  applied  to  the 
numerical  model.  With  IJ ^  =  0, 
then  Vj  =  LJj-V^  =  -V^,  the 
momentum  equation  may  be  w  ritten 
as: 

(14) 

and  now,  taking  advantage  of  our 
mass  conservation  relationship,  we 
have; 


“1 

-V  =  -  V 

w'  T“  0  w 

"0 


We  may  substitute  for  Uq  to  yield; 


=  - 


\S- 


hQ  Hq  + 


hx  \  2 


1/2 


(15) 


(16) 


If  we  consider  the  speed  of  the  perturbation  in  front  of  and  behind  the  shock, 
we  note  that  both  move  toward  the  shock. 

To  demonstrate  this,  we  calculate  the  relative  speeds  Vq  and  Vj.  These  are 
the  speeds  of  fluid  particles  as  perceived  by  an  observer  moving  with  the 
shock.  We  have  already  shown  that  Vy  =  -V^  or 


1  ^0  /. 
T  *  T-  ('■O 


(17) 


The  relative  speed  of  an  upstream  moving  perturbation  Wy  is 

Wi  =  f/j  -  ^  -  1/^  =  Vi  -  (18) 

If  this  value  is  negative,  then  a  perturbation  behind  the  jump  catches  the  shock, 
and  from  Equation  17  we  know  is  greater  in  magnitude  than  VpWj<0. 

In  front  of  the  shock  the  relative  particle  speed  is  Vq. 
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vv»o  = 


(19 


(20) 

(21) 


Again  we  calculate  the  relative  speed  of  a  perturbation,  but  now  in  front  of 
the  shock: 

W'o  =  ^  =  ^  (22) 

Now  if  Wq  is  positive  the  shock  catches  up  with  the  wave  perturbation,  and 
since  Vq  is  clearly  greater  than  this  is  indeed  what  happens.  Therefore 
any  small  perturbations  are  swept  toward,  and  are  engulfed  in  the  shock. 


Shock  relations  in  2-D 

Previous  sections  derive  the 
shock  relations  in  1-D  and  are 
important  for  understanding  behavior 
and  to  produce  test  problems.  Here 
we  extend  these  relations  to  2-D 
(Courant  and  Friedricks  1948).  To 
do  this,  consider  the  region  S2  shown 
in  Figure  4.  It  is  divided  into 
subdomains  and  Qj  by  the  shock 
shown  as  boundary  F^,  which  is 
defined  by  the  coordinate  location 
Xj(t).  The  right  side  boundary  is  F^ 
and  the  left  F|.  The  normal 
direction  is  chosen  as  shown  in 
Figure  4.  Integration  over  the 
subdomains  is  performed  separately; 
and  then  by  letting  the  width  about 
the  shock  go  to  zero,  we  derive  the 
mass  and  momentum  relationship 
across  the  jump  in  the  direction  n. 


Figure  4.  Definition  of  terms  for  2-D 
shock 


Mass  conservation.  For  constant  density  we  have 


£ 

It 


I 


h  da  *  — 
dt 


J  ft  da  =Q 


(23) 
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which  may  be  expanded  as 


/  17  -  f  (V'l  •  -  f  '  1^.(0  •  «J  ^ 

”l  “2  Ti  ^24) 


/*' 

*  ^ 


[^,(0  •  n]  </r  .  V  (V,  • «,)  ^r  =  0 


where, 


V’^=  the  velocity  of  the  left  boundary 

=  the  velocity  of  the  right  boundary 

X^(t)  =  the  velocity  of  the  shock 

h‘  =  the  depth  in  the  limit  as  the  shock  is  approached 
from  subdomain  Qj 

/*■*■  =  the  depth  in  the  limit  as  the  shock  is  approached 
from  subdomain  Q2 

Taking  the  limit  as  ^2  width  we  have 


h'[V-  -  n-h^lV* 


(25) 


where, 

V  =  the  velocity  in  the  limit  as  the  shock  is 
approached  in  subdomain  Qj 

V*"  =  the  velocity  in  the  limit  as  the  shock  is 
approached  in  subdomain  Q2 

For  an  arbitrary  segment  to  preserve  the  equation,  the  integrand  itself 
must  satisfy  the  equation,  therefore 


(26) 
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where 


v;  =  iv  -:f,(0]  •«  (27) 

-  ^:t(01  •«  (28) 

which  states  that  the  relative  mass  flux  jump  across  the  shock  in  the  direction 
n  should  be  zero. 

Momentum  relation.  Again  assuming  constant  density,  the  balance  of 
momentum  and  force  may  be  written  as  (in  the  direction  of  the  normal  to  the 
shock) 


_  r  vh  ‘  n  dQ  *  JL  r 

dt  A  J 


^  Vh  •  n  dQ  ^  -  {  1  ghfn  •  df 


). 

■I. 


(29) 


1  g/j/rt  •  n,  dr 


f  ^(Vh  -n)  dQ  *  f  —  (Vh'  n)  dQ  -  f  (Vj/j,  •  n<)  (V,  •  n)  dT 


*  ^  (y-h-  ^n)  [^,(0  ’n]dr  -  j  n)I2r,(0  •  n]  dT 

"  (30) 


■/, 


{V^r  ■  ”r)(^r  '  n)  dT  ^  -  |  ^  ghfn  •  dT 


I. 


^  J  ghfn  •  /i|  dT 


and  taking  the  limit  as  £2^  and  Q2  shrink  in  width  results  in 
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which  for  an  arbitrary  length  to  preserve  the  equation,  the  integrand  itself 
must  satisfy  the  equation,  therefore; 

-iQ~  •«)v;  *iQ*  •n)v;  *  ^g[  {h  -)2  ]  =  0  (32) 

where, 

Q-  =  V'h* 

Q*  =  V+h"^ 
or 

KQ  (33) 

which  states  that  the  relative  momentum  flux  in  the  direction  n  is  balanced  by 
the  pressure  jump  across  the  shock. 
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2  Numerical  Approach 


The  selection  of  a  numerical  scheme  is  driven  by  two  related  difficulties: 
numerically  modeling  highly  advective  flow  and  the  capturing  of  shocks.  This 
chapter  discusses  the  problem  with  advection  schemes  generally.  It  then 
follows  the  development  of  the  scheme  we  will  use  and  discusses  the 
implications  in  shock  capturing. 


Advection  Dominated  Fiow 

The  problem 

The  quality  of  the  numerical  solution  dq^ends  upon  the  choice  of  the  basis 
(or  interpolation)  function  and  upon  the  test  function.  The  basis  function 
determines  how  the  variable  (or  solution)  is  represented  and  the  test  function 
determines  the  way  in  which  the  differential  equation  is  enforced.  Finite  ele¬ 
ments  are  a  subset  of  the  weighted  residual  method.  Here  one  looks  at  the 
solution  of  a  differential  equation  in  a  weighted  average  sense.  In  the  Galerkin 
approach  the  test  function  is  identical  to  the  basis  function.  This  method  can 
have  difficulty  with  advection-dominated  flow.  The  basic  problem  is  that  the 
form  of  the  test  function  (typically  an  even  or  symmetric  function)  cannot 
detect  the  presence  of  a  node-to-node  oscillation,  since  this  "spurious  solution" 
has  a  spatial  derivative  which  is  an  odd  function  (antisymmetric).  One 
approach  to  resolve  this  problem  is  to  use  a  mixed  interpolation  where,  for  the 
shallow-water  equations,  the  depth  uses  a  lower  order  basis  than  does  the 
velocity  (see,  e.g.,  Platzman  (1^8)  or  Walters  and  Carey  (1983)).  Typically, 
these  are  chosen  as  depth  as  an  elemental  constant  and  velocity  as  linear,  or 
depth  linear  and  velocity  as  a  quadratic.  This  approach  effectively  decouples 
the  depth  from  this  node-to-node  oscillation  but  depends  upon  some  additional 
artificial  viscosity  to  damp  velocity  oscillations  if  the  flow  is  not  highly 
resolved.  Another  approach  is  to  modify  the  test  function  so  that  it  includes 
odd  functions  as  well  as  even  functions  so  that  these  modes  can  be  detected 
and  if  weighted  properly,  eliminated.  Any  approach  in  which  the  test  function 
differs  from  the  basis  function  is  termed  a  Petrov-Galerkin  approach.  In  our 
case  we  choose  the  Lagrange  basis  functions  to  be  C°;  i.e.,  the  functions  arc 
continuous.  Let  us  consider  an  example  to  illustrate  the  problem  with  the 
Galerkin  approach  and  an  approach  to  develop  a  Petrov-Galerkin  test  function. 
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Petrov-Galerkin  formulation 


First  we  will  illustrate  the  problem  that  discrete  formulations  have  with 
advection-dominated  flow.  In  this  regard  the  1-D  linearized  inviscid  Burgers’ 
equation  may  be  written 

C^  =  0  ,  over  domain  L  (34) 


where  the  subscripts  t  and  x  represent  partial  derivatives  with  respect  to  time 
and  space,  respectively,  and 

IJq  ~  the  advection  velocity,  which  here  is  a  constant 

C  =  some  species  concentration 

In  the  discrete  representation  we  shall  approximate  the  solution  as  C*  linear 
Lagrange  basis  functions. 


C{x)  =  2  4)  <jc)  C 
I  ^  ■' 

here  C(x)  is  the  approximate  solution,  and  the  subscript  j  indicates  nodal 
values  and  is  the  Galerkin  test  function  at  node  j. 

Our  numerical  solution  equation,  for  the  steady-state  problem  (Cj=0)  may 
be  written  as  (he  inner  product 


(4),-  ,  Uq  j  4>/  W  Cj)  =0  ,  for  each  i  (35) 


where 


(f(x),  g(x))  -  4  f(x)  g(x)  dx 

and  the  prime  indicates  the  derivative  with  respect  to  x. 

On  a  uniform  grid  the  result  of  this  integration  on  a  typical  patch  is 

^j*l  ~  (36) 


(Note  that  finite  difference  methods  using  central  differences  give  an  identical 
result.) 

In  order  to  demonstrate  that  this  solution  contains  a  spurious  oscillation, 
let’s  write  these  nodal  values  as 
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(37) 


Cj  =  Cpf 


A 

where  C  is  a  constant  determined  by  the  boundary  condition  and  p  is  the 
numerical  root. 

The  roots  of  Equation  36  are 


P 


2 


-1=0 


(38) 


or 


p  =  ±  1 


(39) 


which  makes  the  general  solution 
Cj  =  C  [-Ip' 


(40) 


where  b  is  some  constant. 


The  analytic  solution  corresponds  to  p  =  1.  The  spurious  node  to-node 
oscillation  is  the  root  p  =  -1.  This  root  results  from  a  test  function  which  is 
made  up  solely  of  even  functions;  that  is,  the  test  function,  the  hat  function,  is 
symmetric  about  node  i  (Figure  5).  If  we  consider  the  node-to-node  oscilla¬ 
tion,  its  derivative  is  an  odd  function,  the  inner  product  of  which  with  the  test 
function  is  identically  zero.  This  is  a  solution! 

Now,  if  the  test  function  includes  both  odd  and  even  components,  this 
mode  will  no  longer  be  a  solution.  In  fact,  if  we  weight  the  test  function 
upstream,  these  oscillations  are  damped;  weighting  downstream  amplifies  them. 
A  common  approach  is  to  use  a  test  function,  tp,  weighted  as  follows, 


=  4»,-  +  a  At 


w 

^0 


d<pj 

dx 


(41) 


where  a  is  a  weighting  parameter. 

Here  the  spatial  derivative  supplies  the  odd  component  to  the  test  function. 
The  resulting  discrete  solution  using  this  test  function  is 


\ 


dX 


=  0  ,  for  each  i 


(42) 


from  which  the  numerical  roots  may  be  calculated  by 
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1  C/o  (p^  -  1)  -  a|C/ol  (P  -  if  =0  (43) 

the  roots  of  which  are  then 


p  =  1  ,  - i  (44) 

a  -  _ 

2 

If  a  2  1/2  we  will  have  no  negative  roots  and  therefore  should  not  have  a 
node-to-node  oscillation.  This  spurious  root  that  we  damp  by  increasing  the 
coefflcient  a  is  driven  by  some  abrupt  change,  most  notably  when  some  dis¬ 
continuity  is  required  in  the  equations  due  to  the  imposition  of  boundary 
conditions.  It  is  more  precise  in  a  smooth  region  for  smaller  a. 

The  situation  is  more  complex  for  the  shallow-water  equations,  since  we 
have  a  coupled  set  of  partial  differential  equations.  We  shall  demonstrate  the 
method  used  in  this  model  by  showing  how  it  relates  in  1-D  to  the  decoupled 
linearized  equations  using  the  Riemann  Invariants  as  the  routed  variables. 

The  1-D  shallow-water  equations  in  conservative  form  may  be  written 
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(45) 


dt 


where 


h 


C  = 


P 


P  m  Uh 


F(Q)  = 


p 


If  we  consider  the  linearized  system  with  the  Jacobian  matrix  as  a  constant, 
the  nonconservative  shallow-water  equations  may  be  written  as 

ii.AiC.o 

dl  iX 

where 

0  1 
“  2  2 

-2f/o 

and  the  subscript  0  indicates  a  constant  value. 

We  may  select  the  matrix  P  such  that 

P~^A  P  =  A 

where  A  is  the  matrix  of  eigenvalues  of  A,  and  P  and  P'^  are  composed  of  the 
eigenvectors. 
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A  - 


IXi  0 
0  X2 


f[/o'*’C’o  0 

0  Uq-Cq 


If  we  define  a  new  set  of  variables  (the  Riemann  Invariants)  as 

T  m  PQ 

we  may  write  the  shallow-water  equations  as  two  decoupled  equations 


*  A  il  =  0 


dt 


dX 


(47) 


for  which  it  is  apparent  that  we  can  propose  a  test  function  as 
Jxp  =  I<t>  *  a  AX  IaIA"^ 

dX 


(48) 


which  can  be  returned  to  the  original  system  in  terms  of  the  variable  Q  as 

ftp  =f<t>  *a  AX  P~^  |AjA‘^i»  ^  (49) 

dX 


The  size  and  direction  of  the  added  odd  function  is  then  based  upon  the 
magnitude  and  direction  of  the  characteristics. 

This  particular  test  function  is  weighted  upstream  along  characteristics. 

This  is  a  concept  like  that  developed  in  the  finite  difference  method  of 
Courant,  Isaacson,  and  Rees  (1952)  for  one-sided  differences.  These  ideas 
were  expanded  to  more  general  problems  by  Moretti  (1979)  and  Gabutti  (1983) 
as  split-coefficient  matrix  methods  and  by  the  generalized  flux  vector  splitting 
proposed  by  Steger  and  Warming  (1981).  In  the  finite  elements  community, 
instead  of  one-sided  differences  the  test  function  is  weighted  upstream.  This 
particular  method  in  1-D  is  equivalent  to  the  SUPG  scheme  of  Hughes  and 
Brooks  (1982)  and  similar  to  the  form  proposed  by  Dendy  (1974).  Examples 
of  this  approach  in  the  open-channel  environment  are  for  the  generalized 
shallow-water  equations  in  1-D  in  Berger  and  Winant  (1991)  and  for  2-D  in 
Berger  (1992).  A  1-D  St.  Venant  application  is  given  by  Hicks  and  Steffler 
(1992). 

If  we  analyze  this  approach  on  a  uniform  grid,  we  find  the  following  roots 
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p  =  1 , 


a  X.1  +  i- 
^  2 


a  Xi  -  _ 
'■  2 


(50) 


Again  if  a  ^  1/2,  all  roots  are  non-negative  and  so  node-to-node  oscillations 
are  damped.  In  2-D  we  follow  a  similar  procedure. 

The  particular  approach  to  numerical  simulation  chosen  here  is  a  Petrov- 
Galerkin  finite  element  method  applied  to  the  shallow-water  equations. 


For  the  shallow-water  equations  in  conservative  form  (Equation  1),  the 
Petrov-Galerkin  test  function  ipj  is  defined  as 


iPi-  /  =  4,^  +  a 


y  dx 


*  Ay _ 1  B 

63;  , 


(51) 


where 


a  =  dimensionless  number  between  0  and  0.5 
<t>  =  linear  basis  function 
In  the  manner  of  Katopodes  (1986),  we  choose 

1/2 

1/2 

I  and  q  are  the  local  coordinates  defined  from  -1  to  1. 

A 

To  find  A  consider  the  following; 

A  m  - 

dQ 

P~^  A  P  =  A 

where  A  =  /X  is  the  matrix  of  eigenvalues  of  A  and  P  and  are  made  up  of 
the  right  and  left  eigenvectors. 


Ax  =  2 


Ay  =  2 


dy'^^ 
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A  m  P~^  A  P 


where 


0 


A  = 


0 


0 


0 


0 


^3 


and 


=  U  +  C 

=  u  -  c 

X-,  =  u 

«/ 

c  =  (gh)1^2 

A  similar  operation  may  be  performed  to  define  B. 


Shock  Capturing 

In  the  section,  "Shock  equations,”  in  Chapter  1  we  have  shown  that  unless 
there  is  a  discontinuity  in  depth,  mechanical  energy  will  be  conserved  in  the 
shallow-water  equations  (with  no  friction  or  diffusion).  So  the  obvious  ques¬ 
tion  is  what  happens  in  a  numerical  scheme  in  which  the  depth  is  approxi¬ 
mated  as  C°;  i.e.,  it  is  continuous.  We  are  only  enforcing  mass  and 
momentum,  but  we  are  implicitly  enforcing  energy  conservation.  This  is  the 
result  that  the  GaJcrkin  approach  will  give  using  C°  depth  approximation.  The 
result  is  that  while  mass  and  momentum  con.servation  are  enforced  over  our 
discrete  model,  energy  is  also  conserved  by  including  the  spurious  node-to- 
uodc  mode  we  discussed.  Since  energy  involves  V'^  terms  and  momentum 
only  V,  both  can  be  satisfied  in  a  weighted  average  sense  over  the  region 
included  in  the  test  function.  This  is  due  Ur 
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V^dQ 


(52) 


where  the  lerm  V  means  the  average  value. 

Basically,  energy  is  "hidden"  from  the  numerical  scheme  in  the  shortest 
wavelength  since  the  model  cannot  "see"  this  in  enforcing  momentum  con.scr- 
vation.  So  what  we  need  is  a  scheme  that  damps  this  shortest  wavelength  and 
thus  dissipates  the  energy.  As  we  demonstrated  in  the  previous  section,  this  is 
precisely  what  our  scheme  does.  Therefore,  the  Petrov-Galerkin  scheme  we 
are  using  to  address  advection-dominated  flow  is  a  good  scheme  for  shock 
capturing  as  well.  The  scheme  dissipates  energy  at  the  short  wavelengths. 

We  have  shown  that  when  a  shock  is  encountered,  the  weak  solution  of  the 
shallow-water  equations  must  lose  mechanical  energy.  Some  of  this  energy 
loss  is  analogous  to  a  physical  hydraulic  system  losing  energy  to  heat,  particle 
rotation,  deformation  of  the  bed,  etc;  but  much  of  it  is,  in  fact,  simply  the 
energy  being  transferred  into  vertical  motion.  And  since  vertical  motion  is  not 
included  in  the  shallow-water  equations,  it  is  lost.  This  apparent  energy  loss 
can  be  used  to  our  advantage. 

We  would  like  to  apply  a  high  value  of  a,  say  0.5,  only  in  regions  in  which 
it  is  needed,  since  a  lower  value  is  more  precise.  Therefore,  we  wish  to  con¬ 
struct  a  trigger  mechanism  which  can  detect  shocks  and  increase  a  automati¬ 
cally.  The  method  we  employ  detects  energy  variation  for  each  element  and 
flags  those  elements  which  have  a  high  variation  as  needing  a  larger  value  of 
a  for  shock  capturing.  Note  that  this  would  work  even  in  a  Galerkin  scheme 
since  this  trigger  is  concerned  with  energy  variation  on  an  element  basis  and 
the  Galerkin  method  would  enforce  energy  conservation  over  a  test  function 
(which  includes  several  elements). 

The  shock  capturing  is  implemented  when  Equation  53  is  true 

iJsil  >  Y  (53) 

where 


FTi  -  F 


where  ED^,  the  element  energy  deviation,  is  calculated  by 
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where 


f(E  -  £,)2  dQ 

1/2 

ED.  =  . 

Qi 

j 

£2j  =  element  i 
E  =  mechanical  energy 
aj  =  area  of  element  i 

and  Ej,  the  average  energy  of  element  i,  is  calculated  by 


i  ^ 


and 

E  =  the  average  element  energy  over  the  entire  grid 
S  =  the  standard  deviation  of  all  ED^ 

Through  trial  a  value  of  y  of  1.0  was  chosen. 

An  apparent  limitation  of  this  method  is  that  it  relies  upon  how  the 
elemental  deviation  compares  with  that  of  all  the  other  elements  of  the  grid.  If 
a  problem  contains  no  shocks,  it  would  still  select  the  worst  elements  and  raise 
the  value  of  a.  Conversely,  if  the  domain  contains  numerous  shocks,  it  might 
not  catch  all  of  them.  Perhaps  some  ratio  of  (ED;/^  might  be  meaningful,  and 
should  be  addressed  in  future  studies. 
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3  Testing 


The  testing  of  this  scheme  and  model  behavior  was  undertaken  in  stages. 
These  progress  from  what  is  essentially  a  1-D  test  for  shock  speed  which  can 
be  determined  analytically,  to  a  2-D  dam  break  type  problem  comparison  with 
flume  data,  to  more  general  2-D  geometry  comparison  of  supercritical  transi¬ 
tion  in  a  flume  but  for  steady  state.  This  series  tests  the  model  against  the 
analytic  results  of  the  shallow-water  equations  for  very  limited  geometry,  and 
progresses  to  more  general  geometry  with  the  limitation  of  the  shallow-water 
equations  in  reproducing  actual  flow  problems.  The  applicability  of  the 
shallow- water  equations  to  these  flume  conditions  is  not  so  important  in  this 
study  (since  it  is  interested  in  shock  capturing),  but  is  important  for  model 
application  in  open-channel  hydraulics. 

The  first  test  is  performed  to  determine  the  comparison  of  model  versus 
analytic  shock  speed  in  a  long  straight  flume.  Shock  speed  will  be  poorly 
modeled  if  the  numerical  scheme  is  handled  improperly.  The  analytic  and 
model  tests  are  performed  in  which  the  flow  is  initially  constant  and 
supercritical;  then  the  lower  boundary  is  shut  so  that  a  wall  of  water  is  formed 
that  propagates  upstream.  This  speed  can  be  determined  analytically,  and  a 
comparison  is  made  between  the  analytic  speed  and  the  model  predictions  for  a 
range  of  resolutions  and  time-step  sizes. 

The  second  case  is  a  comparison  to  a  flume  data  set  reported  in  Bell,  Elliot, 
and  Chaudhry  (1992)  which  is  analogous  to  a  dam  break  problem.  Here  the 
shock  is  in  a  horseshoe-shaped  channel  and  the  comparison  is  to  actual  flume 
data.  The  comparisons  are  made  to  the  water  surface  heights  and  timing  of  the 
shock  passage. 

The  final  case  is  a  steady-state  comparison  to  flume  data  reported  in  Ippcn 
and  Dawson  (1951).  Here  a  lateral  transition  under  supercritical  flow  condi¬ 
tions  generates  a  field  of  oblique  jumps.  The  model  comparison  is  made  to 
these  conditions,  which  is  a  more  general  2-D  domain  than  previous  tests. 
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Case  1 :  Analytic  Shock  Speed 


The  shock  speed  for  the  shallow-water  equations  given  simple  1-D 
geometry  can  be  determined  analytically.  These  are  the  Rankine-Hugoniot 
relations  shown  in  Equations  5  and  9.  lliis  provides  a  direct  comparison  with 
the  model  shock  speed  without  relying  upon  hydraulic  flume  data,  for  which 
discrepancy  will  be  due  to  the  hydrostatic  assumption  made  in  the  shallow- 
water  equations.  Instead  we  have  a  direct  way  of  evaluating  the  numerical 
scheme  alone.  As  spatial  and  temporal  resolution  increase,  the  numerical 
shock  speed  should  converge  to  the  analytic  speed.  The  test  consists  of  setting 
a  supercritical  flow  in  a  long  channel,  closing  the  downstream  end,  and 
calculating  the  speed  of  the  jump  that  forms  and  propagates  upstream.  The 
initial  conditions  for  this  test  case  are  shown  in  Table  1.  The  lest  conditions 
are  shown  in  Table  2.  The  term  indicates  the 


Table  1 

initial  Conditions  for  the  Analytic  Shock  Speed  Case 

Condition 

Valiw 

depth 

0.1  m 

U 

(0.3) ’'2 

9 

4  m/s^ 

Table  2 

Test  Conditions  for  the  Analytic  Shock  Speed  Case 

Condition 

Value 

1.0,  1.5 

«.«s 

0.25,0.50 

n 

0.0 

O.O.O.O 

method  applied  to  the  temptoral  derivative,  1.0  is  first-order  backward,  and  1.5 
is  second-order  backward.  The  subscript  s  indicates  the  value  in  the  shock 
vicinity.  The  a  and  are  the  weighting  of  the  Petrov-Galerkin  contribution 
throughout  the  domain  and  in  the  shock  vicinity,  respectively.  With 
Manning’s  n  and  viscosity  of  0.0  there  is  no  dissipation  in  the  shaliow-water 
equations. 

Figures  6-8  and  9-11  show  the  center-line  profile  over  time  of  these  tests 
for  Uj  ~  1.0  and  for  =  0.4  and  0.8  m,  respectively.  These  plots  represent 
the  center-line  depth  profile  over  lime  in  a  perspective  view.  The  vertical  axis 
is  the  flow  depth,  the  horizontal  axis  is  time,  and  the  axis  that  appears  to  be 
into  the  page  is  the  distance  along  the  channel.  From  this  one  can  see  that  as 
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Elevation, 


Figure  7.  Time-history  of  center-line  water  surface  elevation  profiles;  =  1.0,  Ax  =  0.4  m,  At  = 
0.8  sec 
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Elevation, 


Figure  8.  Time-history  of  center-line  water  surface  elevation  profiles;  Oj  =  1.0,  Ax  =  0.4  m,  At  = 
1 .6  sec 


Figure  9.  Time-history  of  center-line  water  surface  elevation  profiles;  Oj  =  1.0,  Ax  =  0.8  m,  At  = 
0.8  sec 
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Figure  10,  Time-history  of  center-line  water  surface  elevation  profiles;  =  1.0,  Ax  =  0.8  m,  At 
=  1.6  sec 


Figure  11.  Time-history  of  (»nter-line  water  surface  elevation  profiles;  =  1.0,  Ax  =  0.8  m,  At 
=  3.2  sec 
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one  moves  over  lime,  the  center-line  profile  shock  moves  upstream.  It  is  apparent  that  as  the 
spatial  and  temporal  resolution  improve,  the  shock  becomes  steeper.  The  shock  is  fairly 
consistently  spread  over  three  or  four  elements;  and  so  as  the  element  size  is  reduced,  the 
resulting  shock  is  steeper.  The  x~t  slope  of  the  shock  indicates  the  shock  speed.  Any  bending 
would  indicate  that  the  speed  changed  over  time,  which  should  not  be  the  case.  The  upper 
elevation  is  precisely  0.2  m,  which  is  correct  There  is  no  overshoot  of  the  jump,  though  there 
is  some  undershoot  when  is  less  than  1.  is  the  product  of  the  analytic  shock  speed  and 
the  ratio  of  time-step  length  to  element  length.  A  value  of  1  indicates  that  the  shock 
should  niove  1  element  length  in  1  time-step. 

Figures  12  and  13  show  the  error  in  calculated  speed  and  the  relative  error  in  calculated 
speed,  re^ectively.  These  are  for  &X  =  0.4,  0.8  and  1.0  m  which  is  reflected  in  the  Grid 
Resolution  Number  defined  as  AY/M.  Here  h  is  the  depth  and  AJt  is  the  analytic  depth 
difference  across  the  shock,  0.1  m.  The  error  was  as  small  as  was  detectable  by  the  technique 
for  measurement  of  ^eed  at  AY  =  0.4  m  so  there  was  no  need  to  go  to  smaller  grid  spacing. 
Values  of  less  than  1  appear  to  lag  the  analytic  shock  and  greater  than  1  leads  the 
analytic  shock.  With  the  largest  the  calculated  shock  speed  is  greater  than  the  analytic  by 
at  most  0.0034  m/sec  which  is  only  0.6  percent  too  fast  As  resolution  is  improved  the 
solution  appears  to  converge  to  the  analytic  ^eed. 

Figures  14-16  and  17-19  are  the  center-line  profile  histories  for  =  1.5  and  for  AY  =  0.4 
and  0.8  m,  respectively.  It  is  apparent  that  the  lower  dissipation  from  this  second-order 
scheme  allows  an  oscillation  which  is  most  notable  upstream  of  the  jump  for  larger  values  of 
Cj.  But  as  Cg  decreases,  there  is  an  undershoot  in  front  of  the  shock.  The  slope  of  the  x-l 
line  along  the  top  of  the  shock  has  a  significant  bend  early  in  the  high  simulations.  The 
speed  is  too  slow  here. 

Now  consider  the  associated  Figures  20  and  21  for  error  in  calculated  shock  speed  and 
relative  error  in  calculated  speed.  The  error  is  actually  worse  than  for  the  first-order  scheme. 
This  is  due  primarily  to  the  slow  ^eed  early  in  the  simulation;  if  this  is  dropped  by  using  only 
the  last  50  seconds  of  simulation,  the  relative  error  is  only  0.6  percent  slower  than  analytic. 
Once  again,  as  the  resolution  improves,  the  solution  converges  to  the  proper  solution. 


Case  2:  Dam  Break 

This  second  case  is  a  comparison  to  hydraulic  flume  results  reported  in  Bell,  Elliot,  and 
Chaudhry  (1992).  A  plan  view  of  the  flume  facility  is  shown  in  Figure  22.  The  flume  was 
constructed  of  Plexiglas  and  simulates  a  dam  break  through  a  horseshoe  bend.  This  is  a  more 
general  comparison  than  Case  1.  Here  the  problem  is  truly  2-D  and  we  now  are  comparing  to 
hydraulic  flume  results,  so  we  must  take  into  consideration  the  limitations  of  the  shallow-water 
equations  themselves.  Initially,  the  reservoir  has  an  elevation  of  0.1898  m  relative  to  the  chan¬ 
nel  bed;  the  channel  itself  is  at  a  depth  (and  elevation)  of  0.0762  m.  The  velocity  is  zero  and 
then  the  dam  is  removed.  The  surge  location  and  height  were  recorded  at  several  stations,  and 
our  model  is  compared  at  three  of  these,  at  stations  4,  6,  and  8.  Station  4  is  6.00  m  from  the 
dam  along  the  channel  center-line  in  the  center  of  the  bend,  station  6  is  7.62  m  from  the  dam 
near  the  conclusion  of  the  bend,  and  station  8  is  9.97  m  from  the  dam  in  a  straight  reach.  The 
model  specified  parameters  are  shown  in  Table  3. 
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Figure  12.  Error  in  model  shock  speed  with  grid  refinement  for  =  1.0 


Figure  13.  Relative  error  in  model  shock  speed  with  grid  refinement  for  a 
1.0 
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Figure  14.  Time-history  of  center-line  water  surface  elevation  profiles;  =  1.5,  Ax  =  0.4  m,  At 
=  0.4  sec 


Figure  15.  Time-history  of  center-line  water  surface  elevation  profiles;  =  1.5,  Ax  =  0.4  m.  At 
=  0.8  sec 
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Figure  16.  Time-history  of  center-line  water  surface  elevation  profiles;  =  1.5,  Ax  =  0.4  m,  At 
a  1 .6  sec 


Figure  17.  Time-history  of  center-line  water  surface  elevation  profiles;  =  1.5,  Ax  =  0.8  m,  At 
=  0.8  sec 
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Figure  18.  Time-history  of  center-line  water  surface  elevation  profiles;  =  1.5,  Ax  =  0.8  m,  At 
=  1.6  sec 


Figure  19.  Time-history  of  center-line  water  surface  elevation  profiles;  cx^  =  1.5,  Ax  =  0.8  m,  At 
=  3.2  sec 
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Table  3 

Test  Conditions  For  The  Dam  Break  Case 

Condition 

Valuo 

“f 

1.0.  1.5 

a,  Oj 

0.25.  0  50 

n 

0.009 

V,  Vj 

0.001 .  0.01  m*/soc 

a 

9.802 

At 

0.05  sec 

The  numerical  grid  is  shown  in  Figure  23,  and  contains  698  elements  and 
811  nodes.  This  grid  was  reached  by  increasing  the  resolution  until  the  results 
no  longer  changed.  The  most  critical  reach  is  in  the  region  of  the  contraction 
near  the  dam  breach.  The  basic  element  length  in  the  channel  is  0.1  m  and 
there  are  five  elements  across  the  channel  width.  For  the  smooth  channel  case, 
Bell,  Elliot,  and  Chaudhry  (1992)  used  a  1-D  calculation  to  estimate  the 
Manning’s  n  to  be  0.016  but  experience  at  the  Waterways  Experiment  Station 
suggests  that  this  value  should  actually  be  0.009,  which  seems  more 
reasonable. 

The  test  results  for  stations  4,  6  and  8  are  shown  in  Figures  24-26.  Here 
the  time-history  of  the  water  elevation  is  shown  for  the  inside  and  outside  of 
the  channel  for  both  the  numerical  model  (at  otj  of  1.0  and  1.5)  and  the  flume. 
The  inside  wall  is  designated  by  squares  and  the  outside  by  diamonds.  Of 
particular  importance  is  the  arrival  time  of  the  shock  front.  At  station  4  the 
numerical  prediction  of  arrival  time  using  ct^  of  1.0  is  about  3.4  sec  which 
appears  to  be  about  0.05  sec  sooner  than  for  the  flume.  This  is  roughly 
1-2  percent  fast  For  of  1.5  the  time  of  arrival  is  3.55  sec  which  is  about 
0.1  sec  late  (3  percent).  At  station  6  both  flume  and  numerical  model  arrival 
times  for  oij  of  1.0  were  about  4.3  sec  and  for  station  8  the  numerical  model  is 
5.6  sec  and  the  flume  is  5.65  to  5.8  sec.  With  Oj  set  at  1.5  the  time  of  arrival 
is  late  by  about  0.2  and  0.15  sec  at  stations  6  and  8,  respectively.  The  flume 
at  stations  6  and  8  has  a  earlier  arrival  time  for  the  outer  wave  compared  to 
the  inner  wave.  The  numerical  model  does  not  show  this.  In  comparing  the 
water  elevations  between  the  flume  and  the  numerical  model,  it  is  apparent  that 
the  flume  results  show  a  more  rapid  rise.  The  numerical  model  is  smeared 
somewhat  in  time  likely  as  a  result  of  the  first-order  temporal  derivative 
calculation  of  otj  of  1.0.  The  numerical  model  with  set  at  1.5  shows  the 
overshoot  that  was  demonstrated  in  Case  1.  This  is  likely  a  numerical  artifact 
and  not  based  upon  physics  even  though  this  looks  much  like  the  flume 
results.  The  surge  elevations  predicted  by  the  numerical  model  are  fairly  close 
if  one  notices  that  the  initial  elevation  of  the  flume  data  is  .supposed  to  be 
0,0762  m  and  it  appears  to  be  recorded  as  much  as  0.015  m  higher  at  some 
gage;.  Since  the  velocity  is  initially  zero  then  all  of  these  readings  should 
have  been  0.0762  m  and  all  should  be  adjusted  to  match  this  initial  elevation. 
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Figure  25.  Flume  and  numerical  model  depth  histories  for  station  6 
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With  this  in  mind,  stations  4  and  8  match  fairly  closely  between  flume  and 
numerical  model.  Station  4  in  the  flume  would  still  have  a  greater  difference 
between  outer  and  inner  wave  than  that  predicted  by  the  model.  The  differ¬ 
ence  might  be  a  manifestation  of  a  three-dimensional  effect  that  the  model 
cannot  mimic.  The  overall  timing  and  height  comparisons  are  good. 

Figure  27  shows  the  spatial  profile  of  the  outer  wall  water  surface  elevation 
of  the  numerical  model  versus  distance  downstream  from  the  dam.  These 
distance  measurements  are  in  terms  of  the  center-line  distance.  The  two  condi¬ 
tions  are  for  of  1.0  and  1.5,  i.e.,  first-  and  second-order  temporal  derivative. 


Outer  Wall  Water  Surface  Elevations 


CbaoDcl  Ceaier  Lioe  Disiance,  m 


Figure  27.  Dam  break  case  water  surface  elevations,  comparison  of 
temporal  representation,  for  time  of  3.5  sec 

The  nodes  are  delineated  by  the  symbols  along  the  lines.  The  overshoot  of 
the  second-order  scheme  and  the  damping  of  the  first-order  is  obvious.  Again, 
it  is  probable  that  the  overshoot  is  a  numerical  artifact  ever  though  this  is 
much  like  what  the  flume  would  show. 


Case  3:  2-D  Lateral  Transition 


TTiis  is  the  most  geometrically  general  case  that  we  test.  The  numerical 
model  is  compared  to  flume  results.  The  flume  data  was  reported  in  Ippen  and 
Dawson  (1951).  The  tests  were  conducted  for  an  approach  Froude  number  of 
4,  upstream  depth  of  0.1  ft,  (0.03048  m)  and  a  total  discharge  of  1.44  ft^/sec 
(0.0408  m^/sec).  The  channel  contracts  from  2  ft  (0.6096  m)  to  1  ft 
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(0.3048  m)  wide  in  a  length  of  4.78  ft  (1.457  m),  i.e.,  an  angle  of  6  deg  on 
each  side. 


The  model  resolution  was  increased  until  we  were  confident  that  the  results 
no  longer  changed  with  greater  resolution.  The  numerical  model  was  set  up 
with  10  evenly  spaced  elements  laterally  across  the  channel  and  24  over  the 
length  of  the  transition.  The  model  limits  were  extended  some  40  ft 
(12.192  m).  The  total  number  of  nodes  was  1661  with  1500  elements.  As  in 
the  flume  test  the  numerical  model  was  set  up  to  provide  a  uniform  depth  of 
0.1  ft  (0.03048  m)  approaching  the  transition.  The  bed  slope  chosen  was 
0.05664.  The  other  parameters  are  shown  in  Table  4. 


Table  4 

Test  Conditions  for  the  Lateral  Transition  Case 

Condidon 

Value 

“/ 

1.0 

a.  ots 

0.10.  0.50 

n 

0.0107 

''■‘'a 

0.005,  0.005  ft^/s  (.0005  m^/s) 

0 

33.208  ft/s*  (9.81 7  m/s^) 

At 

0.005  sec 

Since  the  model  was  run  to  steady-state,  of  1.0  is  appropriate  (time 
accuracy  is  irrelevant  here). 

The  results  from  the  numerical  model  run  and  the  flume  results  are  shown 
in  Figure  28.  The  oblique  shock  forms  along  the  sidewalls  of  the  transition 
and  impinges  on  the  point  in  which  the  converging  channel  goes  back  to  paral¬ 
lel  walls.  This,  by  the  way,  is  the  manner  in  which  one  would  want  to  design 
a  lateral  transition.  The  positive  wave  from  the  beginning  of  the  converging 
walls  will  tend  to  cancel  the  negative  wave  originating  at  the  point  where  the 
walls  change  back  to  parallel.  The  heights  of  the  water  surface  are  indicated 
by  the  contours  in  both  model  and  flume.  The  maximum  and  minimum 
heights  compare  fairly  well.  The  shape  is  good  as  well.  Generally  the  wave 
from  the  shallow-water  equation  will  be  swept  downstream  less  than  that  from 
the  flume  results  since  the  shallow-water  equations  will  transport  all  wave¬ 
lengths  at  the  speed  of  a  long  wave.  Shorter  waves  will  travel  more  slowly 
than  the  shallow-water  equations  predict  The  comparison  is  good,  and  the 
model  demonstrates  that  the  shock  capturing  technique  functions  well  in  a 
general  2-D  setting. 
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DISTANCE  FROM  CENTER  LINE,  FT  DISTANCE  FROM  CENTER  LINE,  FT 


Figure  28.  Comparison  of  flume  and  numerical  model  water  surface  elevations  for  the  super¬ 
critical  transition  case,  straight-wall  contraction  F  =  4.0.  To  convert  feet  to 
meters,  multiply  by  0.3048 
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Discussion 

Now  let’s  study  the  behavior  of  the  1-D  linearized  shallow-water  equation 
analytically  and  numerically.  This  could  lead  to  a  concqjtual  appreciation  of 
the  behavior  we  have  observed  in  the  testing  section  of  the  report.  We  shall 
follow  a  Fourier  analysis  of  the  wave  components;  for  examples,  see 
Leendertse  (1967)  or  Froehlich  (1985).  First  let's  consider  the  nondimension- 
alized  shallow-water  equations 


(54) 


where,  the  subscript  *  indicates  nondimensional  quantities  and  o  as  a  subscript 
indicates  a  constant,  and 


\h/h^ 


These  equations  can  be  diagonalized  by  defining  a  new  variable 


T  =  P~^Q^ 


(55) 


such  that  P^A^Po  = 


where 


A  o 


0  1 

'1-^1 


Po 


0 

O' 

0 

Po-\ 

0  X2 

^2  -1 
-1 
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Ag  is  the  diagonal  matrix  of  eigenvalues  and  and  P'^  are  composed  of  the 
eigenvectors  (and  are  arbitrary).  With  the  substitution  of  Equation  55  into  54 
and  multiplication  by  we  retrieve  the  diagonalized  shallow-water  equations 
in  terms  of  the  Riemann  Invariants 


dT 

K 


(56) 


Now  if  we  consider  solutions  in  terms  of 

T{x,t)  =  fe'<^e^  (57) 

A 

where  7  is  a  constant  and  AT  is  the  wave  number,  we  arrive  at  the  solution 

r(x,r)  =  f  e‘(^^^)  (58) 


where 

CO  =  Kk,  the  wave  frequency 
Y  =  -ICO 

With  this  solution  we  shall  now  compare  the  behavior  of  the  model  to  tha^  of 
the  analytic  solution. 

The  test  function  for  Equation  54  in  HIVEL2D  would  be 

nii  /  =  ♦,/*  oA.  p;'  |A„i  a;’  ^  m 

Now,  since  7  is  a  linear  combination  of  the  variables  h*  and  P,  we  can  con¬ 
vert  this  to  the  diagonal  system  as  well,  so  that  the  equivalent  test  function  is 

+  oAx  lA^lA'^  (60) 


implying  this  test  function  to  the  discretized  differential  equation  and 
substituting 


m  =  Y:<l,j(x)rj 
j 


and 
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y.  «  -  -intdii)  gij/SxK 

where  the  superscript  n  indicates  the  time-step  and  the  subscript  j  is  the  spatial 
node  location. 

We  now  present  the  results  of  this  analysis  for  a  =  1/2  and  for  the  temporal 
derivative  parameter  of  1.0  and  1.5.  We  shall  compare  the  relative  ampli¬ 
tude  and  relative  speed  for  a  single  time-step.  The  parameter  for  relative  ^eed 
is  given  by 


relative  speed  - 


N 


(61) 


where 


I 


i 


N  -  elements  per  wavelength 

XAt* 

C  =  Courant  number  ■  _ 

Ax, 

X  =  wave  speed,  either  or  X2 

For  ttj  =  f ,  which  is  first-order  backward  difference  in  time,  the  relative 
amplitude  is  shown  in  Figure  29  and  the  relative  wave  speed  is  shown  in  Fig¬ 
ure  30.  This  is  plotted  versus  the  number  of  elements  per  wavelength  N  and 
the  Courant  number  C.  Also  remember  that  these  comparisons  apply  for  either 
characteristic  (X^  or  X^),  even  for  subcritical  conditions  in  which  X2  is 
negative.  In  these  figures  the  Courant  number  varies  from  0.5  to  2.0  and  the 
elements  per  wavelength  from  2  to  10. 

The  amplitude  portrait  shows  substantial  damping  for  larger  C  and  for  the 
shorter  wavelengths  (or  alternatively  the  poorer  resolution).  The  large  damp¬ 
ing  at  a  wavelength  of  2Ajt  is  important,  as  this  is  the  mechanism  that  provides 
the  energy  dissipation  to  capture  shocks.  Now  consider  the  phase  portrait,  or 
in  this  case  the  relative  speed  portrait.  Over  the  conditions  shown,  the  numeri¬ 
cal  speed  is  less  than  the  analytic  speed  throughout  For  larger  C  the  relative 
speed  is  somewhat  lower  (worse).  For  A/  =  2  the  speed  is  0,  so  that  undamped 
oscillation  could  remain  at  steady  state. 
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In  comparison  to  the  results  we  have  shown  in  Figures  6-11  for  Case  1, 
analytic  shock  case,  we  must  remember  that  is  the  Courant  number  based 
on  shock  speed,  whereas  C  is  based  on  the  perturbation  wave  speed.  If  we 
consider  a  wave  moving  upstream  just  behind  the  shock,  since  short  wave¬ 
lengths  move  too  slowly,  the  disturbance  of  the  shock  produces  waves  of  these 
length  which  fall  behind  the  shock  rather  than  remaining  within.  As  the  time- 
step  is  reduced  (C^  gets  smaller)  the  relative  speed  is  better  for  the  moderate 
wavelengths  and  so  the  shock  front  becomes  sharper. 

At  a  point  near  the  shock  front  we  note  that  generally  we  get  a  sharp  front 
with  no  undershoot  until  we  reach  the  smallest  time-step.  Again  if  we  are 
within  the  shock  at  a  depth  where  there  is  an  upstream  propagating  wave 
(subcritical),  is  there  a  Courant  number  C  that  has  a  relative  speed  greater  than 
analytic.  This  would  be  the  only  way  in  which  an  undershoot  could  appear. 
Figure  31  extends  the  relative  wave  speed  portrait  below  C  =  0.5.  From  this 
figure  it  is  apparent  for  small  values  of  C  that  the  numerical  wave  speed  is 
greater  than  analytic  so  th"'*  it  is  possible  to  develop  an  undershoot  in  front  of 
the  jump. 

For  =  1.5  we  have  a  second-order  temporal  derivative  which  has  relative 
amplitude  and  relative  speed  portraits  shown  in  Figures  32  and  33,  respec¬ 
tively.  The  degree  of  damping  is  much  less  than  for  the  first-order  case.  The 
relative  speed  is  better  but  not  so  dramatic  as  the  improvement  in  amplitude. 
An  interesting  point  is  that  the  relative  speed  for  ^  =  2  is  nonzero  for  lower  C 
values.  This  implies  that  a  spurious  mode  should  not  reside  in  the  grid  at 
steady  state.  In  Figure  34,  we  show  the  relative  ^eed  portrait  extended  below 
C  values  of  0.5.  As  with  =  1,  for  very  low  C  the  numerical  relative  speed 
is  greater  than  the  analytic.  Therefore,  we  would  expect  to  have  an  undershoot 
for  small  time-steps.  It  should  become  more  pronounced  and  longer  as  the 
time-step  is  reduced  further.  Since  we  generally  have  a  relative  speed  lower 
than  analytic,  we  expect  an  overshoot  behind  the  jump  which  becomes  longer 
as  the  time-step  is  increased.  Refening  to  Figures  14-19  of  case  1,  this  is 
precisely  what  we  note.  Also,  for  smaller  time-steps  there  is  some  undershoot 
as  well.  These  same  features  are  notable  in  the  second  test  case,  the  dam 
break  test  case. 

For  the  sake  of  completeness  the  relative  amplitude  and  speed  portraits  are 
included  for  a  =  0  and  0.25  at  of  1.0  and  1.5  in  Figures  35-42.  The  condi¬ 
tion  a  =  0  is,  in  fact,  the  Galerkin  case  since  the  Petrov-Galerkin  contribution 
is  included  through  a.  The  Galerkin  approach  is  shown  to  contain  a  steady- 
state  spurious  mode  due  to  the  speed  of  zero  for  N  =  2.  Furthermore,  this 
mode  is  undamped.  The  case  of  a  =  0.25  shows  that  the  relative  speed 
portraits  change  very  little  from  a  =  0.5  but  the  amplitude  damping  is 
improved. 

The  obvious  conclusions  that  can  be  drawn  from  this  discussion  is  that  for 
an  unsteady  run  either  use  =  1.5  or  take  smaller  time-steps  with  =  1.0. 

An  improvement  in  spatial  resolution  dramatically  improves  the  solution. 
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Figure  31.  Relative  speed  versus  C  and  resolution  tor  =  1.0  and  a  =  0.5,  for 
small  values  of  C 
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Figure  33.  Relative  speed  versus  C  and  resolution  for  =  1.5  and  a  =  0.5 


Figure  34.  Relative  speed  versus  C  and  resolution  for  =  1.5  and  a  =  0.5,  for 
small  values  of  C 
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Figure  35.  Relative  amplitude  versus  C  and  resolution  for  =  1.0  and  a  =  0 


50 


Chapter  3  Testing 


Figure  37.  Relative  amplitude  versus  C  and  resolution  for  a.  -  1.0  and 
a  =  0.25 


Figure  38.  Relative  speed  versus  C  and  resolution  for  =  1.0  and  a  =  0.25 
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Figure  39.  Relative  amplitude  versus  C  and  resolution  for  a,  =  1 .5  and  a  =  0 


Figure  40.  Relative  speed  versus  C  and  resolution  for  =  1.5  and  a  =  0 
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Figure  41 .  Relative  amplitude  versus  C  and  resolution  for  =  1.5  and 
a  =  0.25 


Figure  42.  Relative  speed  versus  C  and  resolution  for  =  1.5  and  a  =  0.25 


ChaiiUtr  3  Tasting 


4  Conclusions 


In  this  report  an  algorithm  is  developed  to  address  the  numerical  difficulties 
in  modeling  surges  and  jumps  in  a  computational  hydraulics  model.  The 
model  itself  is  a  finite  element  computer  code  representing  the  2-D  shallow 
water  equations. 

The  technique  developed  to  address  the  case  of  advection-dominated  flow  is 
a  dissipative  technique  that  serves  well  for  the  capturing  of  shocks.  The 
dissipative  mechanism  is  large  for  short  wavelengths,  thus  enforcing  energy 
loss  through  the  hydraulic  jump,  unlike  a  nondissipative  technique  used  on  C” 
representation  of  depth,  which  will  implicitly  enforce  energy  conservation, 
dictated  by  the  shallow-water  equations,  through  a  2Ax  oscillation. 

The  test  cases  demonstrate  that  the  resulting  model  converges  to  the  correct 
heights  and  shock  speeds  with  increasing  resolution.  Furthermore,  general  2-D 
cases  of  lateral  transition  in  supercritical  flow  showed  the  model  to  compare 
quite  well  in  reproducing  the  oblique  shock  pattern. 

The  trigger  mechanism,  based  upon  energy  variation,  appears  to  detect  the 
jump  quite  well.  The  Petrov-Galerkin  technique  shown  is  an  intuitive  method 
relying  upon  characteristic  speeds  and  directions  and  produces  a  2-D  model 
which  is  adequate  to  address  hydraulic  problems  involving  jumps  and  oblique 
shocks. 

The  resulting  improved  numerical  model  will  have  application  in  supercriti¬ 
cal  as  well  as  subcritical  channels,  and  transitions  between  regimes.  The 
model  can  determine  the  water  surface  heights  along  channels  and  around 
bridges,  confluences,  and  bends  for  a  variety  of  numerically  challenging  events 
such  as  hydraulic  jumps,  hydropower  surges,  and  dam  breaks.  Furthermore, 
the  basic  concepts  developed  are  applicable  to  models  of  aerodynamic  flow 
fields,  providing  enhanced  stability  in  calculation  of  shocks  on  engine  or  heli¬ 
copter  rotors,  for  example,  as  well  as  on  high-speed  aircraft. 
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